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Abstract 

Apoptosis is one of the most basic biological processes. In apoptosis, tens of species are involved 
in many biochemical reactions with times scales of widely differing orders of magnitude. By the 
law of mass action, the process is mathematically described with a large and stiff system of ODEs 
(ordinary differential equations). The goal of this work is to simplify such systems of ODEs with 
the PEA (partial equilibrium approximation) method. In doing so, we propose a general framework 
of the PEA method together with some conditions, under which the PEA method can be justified 
rigorously. The main condition is the principle of detailed balance for fast reactions as a whole. 
With the justified method as a tool, we made many attempts via numerical tests to simplify the 
Fas-signaling pathway model due to Hua et al. (2005) and found that nine of reactions therein can 
be well regarded as relatively fast. This paper reports our simplification of Hua at el.'s model with 
the PEA method based on the fastness of the nine reactions, together with numerical results which 
confirm the reliability of our simplified model. 

Keywords: Partial equilibrium approximation, apoptosis, biochemical reactions, the principle of 
detailed balance, sensitivity analysis 

1 Introduction 

Apoptosis is one of the most basic biological phenomena. It is a cellular suicide route that allows for 
the selective removal of superfluous and potentially dangerous cells. This genetically controlled process 
ensures normal embryonic development, tissue homeostasis and normal immune-system function in mul- 
ticellular organisms. On the other hand, defects in apoptosis may cause serious diseases such as cancer, 
autoimmunity, and neurodegeneration [TJ [SJ El S] . For these reasons, understanding the mechanism of 
apoptosis is of fundamental importance. 

The apoptotic process involves tens of biological molecules (species) , which react within tens of bio- 
chemical reactions with time scales of widely differing orders of magnitude. When the law of mass action 
[5] is employed, it is described mathematically by a simultaneous system of tens of ordinary differential 
equations (ODEs). Such a large scale and stiff system of ODEs can hardly help us to understand the 
mechanism of the apoptosis. 
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The goal of this work is to derive mathematically reliable simplifications of the large apoptosis system 
proposed by Hua et al. in [6 for human Jurkat T cells. Two widely used methods for simplifying chemical 
kinetics are the Quasi Steady-State Approximation (QSSA) [5], also called the Bodenstein method, 
and Partial Equilibrium Approximation (PEA) [9l [TUJ [11] . The former assumes that the concentrations 
of transient intermediate species reach steady states and thereby the rate equations for the intermediate 
species are replaced with algebraic relations. On the other hand, the PEA simply takes the fast reactions 
in equilibrium. In this manner, the stiffness is removed and some algebraic constraints are obtained. 
For both methods, the algebraic relations can be used to reduce the number of the rate equations and 
consequently the chemical kinetics is simplified. An unexpected benefit of such simplifications is that 
less parameters are needed for the simplified models than for the original ones. This is good because the 
parameters are often not reliably known (see [9]). 

The QSSA, PEA and their combinations have been extensively used to simplify chemical kinetics 
mechanisms for many years, with great success [TH [T21 HH [HI EH HZJ [TB] . They were also used by Okazaki 
et al. in [T5] to simplify the large apoptosis system in [BJ (see comments below). However, these methods 
seem to lack a systematically mathematical justification. Recently, we pointed out that the PEA method 
can be rigorously justified for reversible reactions obeying the principle of detailed balance [201 [H], by 
using the singular perturbation theory of initial- value problems for ODEs |22) . Thus, our simplification 
will base on this justified PEA method and therefore is reliable. 

As commented in [25] , Okazaki et al.'s simplification seems baseless. In fact, when applying the 
QSSA method for the intermediate species Casp82:Casp3, Okazaki et al. assumed that the concentration 
sum of the intermediate species and Casp82 (an activated initiator caspase) was conserved and obtained 
the Michaelis-Menten equation for the product (see Appendix A of [H]). The latter is only true if the 
activated initiator caspase does not participate in other reactions. However, this is not the case here 
because the activated caspase is simultaneously, instead of consecutively, involved in other reactions. 
Similar conservation assumptions were used for several steps. This is why we question the simplified 
model in [19] . although it is successful in some sense. 

This paper is a continuation of our previous work [23] . In |23| . we showed that two molecules (Smac 
and XIAP) involved in the apoptotic process, neglected in [T2], are not negligible in general. Then we 
applied the justified PEA method to obtain a very preliminary simplified model by assuming only six 
reversible reactions to be fast. In the present paper, we will use the PEA method to further simplify 
the large apoptosis system [BJ. To do this, we firstly verify the principle of detailed balance for more 
fast reactions as a whole. Having such a verification, the singular perturbation theory of initial-value 
problems for ODEs can be employed to derive our simplified models. Then we use numerical simulations 
to compare the new models with the original one and Okazaki et al.'s simplified model from various 
aspects, including accuracy, sensitivity and the M-D transition behavior [19] . Moreover, we introduce 
a new quantity to evaluate our simplifications. All numerical results confirm the reliability of both our 
simplified models and the PEA method. 

Let us remark that, thanks to its reliability, the justified PEA could be used as a tool to determine 
whether or not a reversible reaction is relatively fast. To this end, one could numerically compare the 
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solution of partial equilibrium computation with that of the fully non-equilibrium computation. If the 
two solutions are close to each other, the reaction can be claimed to be relatively fast. 

The paper is organized as follows. In Section 2 we present a general framework of the PEA method, 
together with some conditions under which the method can be justified rigorously. The apoptosis process 
is introduced in Section 3. The PEA method is used in Section 4 to simplify the apoptosis system by 
checking the principle of detailed balance for nine reversible reactions as a whole. Numerical simulations 
are reported in Section 5. Finally, the main results of this paper are summarized in Section 5. 



2 The PEA Method 

In this section we first present a general framework of the PEA method, together with some conditions 
under which the method can be justified rigorously. Then we take the simplest system for enzyme 
inhibition as an example to show how to use our PEA method. 

2.1 A general framework of the PEA method 

Consider a system with N chemical species Ci(i = 1, 2, • • ■ , N) participating in M reactions 

aid + a\C 2 + ■■■ + a p N C N ; = > 6?d + b%C 2 + --- + b p N C N (2.1) 

k—p 

for p = 1, 2, • • • , M . Here the non- negative integers a\ and b\ are the stoichiometric coefficients of the 
z tft, -species in the p t,l -reaction, and fc +p and fc_ p are the respective forward and backward rate constants 
of the p* ft -reaction. The reversibility means that both k +p and k- p are positive. 

Denote by u% — [Ci](t) the concentration of the z tft, -species C\ at time t. According to the law of 
mass action [5], the evolution equation for u, is 

— -^K-<H (2-2) 



with 



dt ^ 

P =i 



_ U °1 a 2 a N _ J. b l b 2 b N 



being the reaction rate of the p-th reaction. 

Suppose the first M'(< M) reactions in (|2.1[) are much faster than others. Then the kinetic equations 
in (|2.2p can be rewritten in the vectorial form: 

% = -Qi(y) + Q 2 (y,z), 

dt s (23) 

Here V is an iV'-vector consisting of those Ui so that the i^-species participates in the fast reactions, Z 
consists of the rest itj, e is a small positive parameter characterizing the fastness, and Qx(V), Q 2 (V, Z) 
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and Qs(V,Z) stand for the corresponding reaction rates. The introduction of e is to make Qi(V) have 
the same order of magnitude as Q 2 {V, Z) and Qz(V, Z). A special case is that Z is void and V contains 
all Ui. 

Assume that there is a steady state u* — (u*,U2, ■ ■ -u^) satisfying u* > for all i and a zero net 
flux condition for each fast reaction: 

k +p (ul)<(u* 2 ) a i ■ ■ ■ (u* N )< - fc_ p «) 6 ?K) b ? ■ ■ • (u* N fi = 0, Vp = 1,2,-.. ,M'. 

This is just the principle of detailed balance 20 for the partial system consisting of the fast reactions 
merely. Clearly, it can only be true when both k +p and fc_ p are positive, that is, reversible reactions. 

Under this assumption, we know from |21) that there is a strictly convex function r/ = rj{V) so that 
Qi(V) can be written as 

Qi(V) = S(V) VV (V) 

for V with strictly positive components. Here S(V) is a symmetric matrix with null-space independent 
of V and r]v(V) is the gradient of r)(V). Moreover, the singular perturbation theory [22] for initial- value 
problems of ODEs can be applied to the stiff system in (|2.3D . In particular, the solutions to initial- value 
problems of (|2.3I) converge uniformly to those of a corresponding reduced system, as e goes to zero, in 
any bounded time interval away from zero. 

In order to derive the reduced system, we notice the ^-independence of the null-space and denote 
by II the constant matrix whose rows span the left null-space of S(V). Without loss of generality, we 
assume that II is of the form 

n = (e,/) 

with I the unit matrix of proper order. Accordingly, we introduce the partition 

v=( X ), cM^f^'V Q, { v,z)=[^ Y - z) ). 

This partition ensures that X can be uniquely and globally obtained by solving Qi(X, Y — QX) = (see 
[2T] if necessary). 
Define 

Y = Y + ex. 

The kinetic equations in (|2.3|) can be rewritten as 

= -Q 1 (X,Y-QX) + Q 1 (X,Y-QX,Z), 
at e 

dY 

— = Q 2 (X, Y - QX, Z) + 6Qi(X, Y - QX, Z), (2.4) 
at 

^L = q 3 (x,y -ex, z). 
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As e goes to zero, the reduced system for (|2.4I) is 



Qx(X,Y-QX) = 0, 

— = Q 2 (X, Y - QX, Z) + 6Qi(X, Y - 9A, Z), (2.5) 

— = Q 3 (x,Y-ex, Z). 

From Qi(X, Y — QX) = we solve X in terms of Y, say X = $(Y). Substituting this expression into 
the second and third equations in (|2.5[) . we obtain 

' dY - ~ ~ - 

— = Q 2 ($(Y),Y - e$(Y), Z) + QQxiHY), Y - Q$(Y),Z), 

This is our simplified model by using the PEA method. Note that the original variables X and Y are 
recovered from 

x = $(Y), y = r-e$(T). 

In case X can be solved from Qi(X,Y) — in terms of Y, say A = ^(Y), we recall the fact that X can 
be obtained by solving Qi(X, Y — QX) = and may well assume that both the Jacobian matrices Q\x (of 
Ql(X,Y) with respect to A) and [Qix-QiyQ] are invertible. Then [I-Q~ X Q 1Y Q] = Qi X [Qix -Qiy@] 
is invertible. It is an easy exercise to show that the invertibility of [/ — Qi X QiyQ] is equivalent to that 
of [I - QQixQiy}- On the other hand, we deduce from Qi(*(Y),Y) = that Q 1X ^(Y) Y + Qiy = 
and thereby *(Y)y = -Qi X Qiy- Now we compute from Y = Y + 0*(Y) that 

Thus we gain equations for Y: 
dY 

— = (I + ey(Y) Y )- 1 (Q 2 (y(Y),Y,Z) + QQ 1 (y(Y),Y,Z)). 

Consequently, the reduced system can be written as 
dY 

— = (i + e*(r) 1 ; 1 )(Q 2 (vi/(r), y, z) + qq 1 (^(y), y, z)), 

^- = Q^(Y),Y,Z) 



(2.6) 



together with the algebraic relation X = \&(Y). 



2.2 A simple example 

In order to elucidate how to use the PEA method, we consider the simplest system for enzyme inhibition 
[S]. This system is demonstrated graphically in Fig. [T]and reveals the competitively inhibitory mechanism, 
where the enzyme reaction is stopped when the inhibitor is bound to the active site of the enzyme. In 
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E + I^C 2 . 



Figure 1: The simplest mechanism for enzyme inhibition 



Fig. [TJ the symbols E, S, I, P, C\ and C 2 stand for the enzyme, the substrate, the inhibitor, the product 
and two complexes, respectively. k\ and fc_i are the respective kinetic rate constants of the forward and 
backward reaction for substrate binding, while and £;_3 are those of the inhibitor binding reaction. k% 
is the rate constant of the substrate conversion reaction. 

According to law of of mass action, the corresponding kinetic equations read as 

d[Ci] 

- — vi - v 2 



dt 
d[C 2 ] 

dt 
d[E] 



v 3 

—vi - v 3 + v 2 



-dT = - Vl 

m 

dt 

d[P] 



dt ~ V3 



dt ^ 
with 

v 1 = k 1 [E\[S\-k- 1 [Ci], 
v 2 = k 2 [Ci], 

v 3 = k 3 [E][I}-k- 3 [C 2 ]. 

Classically, the reactions for enzyme to bind to substrates and inhibitors are regarded as fast and 
reversible. Thus we rewrite (|2.7|l as 



~V\ — v 2 

1)3 



djd] = 1 

dt e 

d[C 2 ] = 1 

dt s 

d[E] _ 1 A 1 

dt e 

d[S] 1 A 

dt e 

d[I] 1 . 
—rr = — "3 
dt e 

d\P] 



-Vl V3 + v 2 



(:» 



where v\ = ewi, W3 = ev^ 1 and v\ and V3 have the same order of magnitude as v 2 . Thanks to the 
reversibility, it is obvious that there are positive numbers [E]*, [5]*, [Ci]*, [/]* and [C2]* such that 



i>i = h[E}*[S}* - fc_i[Ci]* =0, t> 3 = fc 3 [£]W - fc_ 3 [C 2 



0. 



Thus the PEA method above can be well applied to the stiff system 
The reduced system for (|2.8|) can be derived as follows. Set 



X= ([Ci],[C a ]) , 
y= ([£7],[5],[i]) T 1 
Z= [P]. 



the stiff system (|2.8p can be rewritten as 



dX_l I vi 
dt e 



V-2 



!'3 







d(Y + QX) 



dt 



t'2 



dt 



(° ) 

-1 



V-2 



with 
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1 

V o i; 



From «i = U3 = 0, we get 



[C 1 ]=X 1 [S][5], [C 2 ]=K 3 [E][I] 



(2.9) 



(2.10) 



with = kj/k-j for j = 1,3. Substituting these into the last two equations in (|2.9[) , we obtain the 
following ODEs 



d([E] + [d] + [C 2 ]) 
dt 

d([S] + [C\]) 



= 0, 



dt 

d([I} + [C 2 ]) 



dt 



= -h i Kx\E\\S\, 
= 0, 



d\P\ 
dt 



= k 2 K 1 [E}[S}. 



This, together with (|2.10l) . is the simplified model by using the justified PEA method. This model can 
be further simplified by using the conservation laws indicated in the first and third equations. We omit 
it here and leave it to the interested reader. 
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3 Apoptosis Systems 

Here we introduce the large apoptosis system proposed by Hua et al. in [5J for human Jurkat T cells. 
To begin with, we recall that there are at least two pathways to trigger apoptosis — intrinsic (mito- 
chondrial) and extrinsic (death receptor) signalling pathways. Both induce death-associated proteolytic 
and/or nucleolytic activities. The intrinsic pathway is initiated when the cell is severely damaged 
or stressed, while the extrinsic one is activated when extracellular death ligands are bound by their 
cognate membrane-associated death receptors such as TNF-Rl(DRl,p55), Fas(DR2,CD95), DR3(APO- 
3,TRAMP), DR4(AP0-2,TRAIL-R1) and DR5(TRICK2,TRAOL-R2) [IJ E2 ESI ESI EI] • 

The Fas-induced signaling pathway is among the best understood and can be schematically shown in 
Fig. El It begins with the binding of Fas ligands (FasL), Fas and FADD (Fas- associated death domain) 
to form the complex DISC (death-inducing signaling complex). The latter can recruits initiator caspases 
such as caspase-8 (Casp8) molecules to cleave and activate them. The activated initiator caspase (CaspS^) 
can cleaves and activates the executor caspase-3 (Casp3) to form Casp3* directly. The amount of Casp3* 
is the indicator of apoptosis. This way to activate Casp3 is called D-channel. In addition, Casp3 can 
also be activated in a so-called M-channel. In this channel, Casp8*. cleaves Bid to generate truncated 
(t)Bid. The tBid then binds to two molecules of Bax to form a complex tBid:Bax2, which will induce 
the release of Cyto.c and Smac from the mitochondria. The released Cyto.c* will combine an adaptor 
protein Apaf-1, ATP and caspase-9 to form apoptosome and thereby activate caspase-9. The activated 
caspase-9 (Casp9*) cleaves and activates Casp3. On the other hand, the M-channel can be blocked by 
XIAP (X-linked inhibitor of apoptosis protein) and Bcl2 through their bindings to the released Smac*, 
Casp9, Casp3*, Bax and tBid. 

The Fas-signaling pathway model proposed by Hua et at. [5] consists of biochemical reactions (Hl)- 
(H25) given in Table 1. From this table we see that the process activating the initiator caspase-8 (Casp8) 
consists of the reactions from (HI) to (H6), which is initiated by FasL. The activated Casp8*> enzymatically 
cleaves caspase-3(Casp3) to produce activated executor Casp3* ((H7) and (H8)) and Bid to generate tBid 
((H9) and (H10)) simultaneously. Then the tBid associates with two Bax to form tBid:Bax2 through 
(Hll) and (H12), which induces the release of Cyto.c and Smac from mitochondrial to cytosol ((H15) 
and (H13)). The released Cyto.c (Cyto.c*) combines Apaf (Apaf-1) and ATP to form an apoptosome 
(Cyto.c* :Apaf: ATP) in (H16), which recruits two caspase-9(Casp9) and generates the activated caspase-9 
(Casp9*) through (H17) to (H19). The activated Casp9* can also enzymatically cleaves and activates 
caspase-3 ((H20) and (H21)). On the other hand, the roles of Casp9, Casp3*, tBid and Bax can be 
inhibited by binding to XIAP((H22) and (H23)) and Bcl 2 ((H24) and (H25)), while the released Smac 
(Smac*) can suppress the function of XIAP (H14). Table 1 also contains all the forward/backward rate 
constants k±i(i = 1,2, • • • , 25). 

Observe that not every reaction in Table 1 is reversible and the reactions activating Casp8 are 
independent of the rest. As in [19 , we call the downstream process, consisting of reactions from (H7) 
to (H25), as the intracellular-signaling subsystem (ISS). Moreover, we follow [19] and assume that the 
concentration of ATP is a fixed constant. Thus, there are 28 species and 19 biochemical reactions involved 
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Table 1: The Fas-signaling pathway model due to Hua et at. (2005) 



(HI) 
(H2)° 
(H3) b 
(H4) b 

(H5) c 

(H6) 

(H7) 

(H8) 
(H9) 

(H10) 
(Hll) 

(H12) 

(H13) 
(H14) 

(HIS) 
(H16) 
(H17) 
(HIS) 

(H19) 
(H20) 

(H21) 
(H22) 

(H23) 

(H24) 

(H25) 



FasL 4- Fa 



, FasC 



9.09 X 10" 



FasC : FAD Dp : Casp8 q : F LI P r 4- FADD 
FasC : F AD Dp : CaspSg : F I L P r + Casp8* 
FasC : FADDp : Casp8 q : F L I P r 4- FILP^ 



k -H2 
k H3 

k -H3 
k H4 



■-Fas : FADD p+1 : Coj P 8, : F I L P r 
t Fas : FADD p : Casp8 q+1 : F I L P r 
-Fas : FADDp : Casp8 q : FILP r + 1 



k H5 



FasC : FADDp : Casp8 q : FL 
k H6 

Casp&q : p41 t-C asp8^ 

OaspH + Casp3 ; • CaspS : Casp3 



*Casp8* 2 ■■ P41 + FasC : FADDp : Casp8 q _ 1 : 



* : Ca.p3- 



k -H7 
k H8 



*-Casp8* + Casp3* 



k H9 

Casp8* + Bid* 1 Casp8* : Bid 



K H10 

Casp8* : Bid ►CaspS* + tBid 



-Hll 



tBid : Bax 4- Bax ^ 



V H12 



i tBid : Bax 2 
k -H12 
k H13 

: 4- tBid : Bax 2 *-S- 



■* 4- tBid : Bax 2 
XI AP 



k H15 

Cyto.c + tBid : Bax 2 •■Cyto.c* + tBid : Bax 2 

k H16 

Cyto.c + Apaf + ATP == * Cyto.c : Apaf : AT P 

k -H16 

k H17 

Cyto.c* : Apaf : ATP + Casp9* 1 Cyto.c* : Apaf : ATP : Casp9 



Cyto.c* : Apaf : ATP : Casp9 + Casp9 =; 



-H18 



Cyto.c* : Apaf : ATP : Casp9 2 



Cyto.c* : Apaf : ATP : 
k H20 



k H19 



*Cyto.c* : Apaf : ATP : Casp9 + Casp9* 



Casp9* + Casp3 5 



•-Casp9* : Casp3 



k -H20 
k H21 

Casp9 : Casp3 *Casp9 + Casp3 

Casp9 + XIAP . kH22 . Casp9 : XI AP 



Casp3* + XIAP 



-H22 
k H22 



k -H22 
k H24 

Bcl 2 + Bax =± Bcl 2 : Bax 

k -H2i 
k H2b 

Bcl 2 + tBid =1 Bcl 2 : tBid 



Casp3* : XIAP 



n-4„ 



3.50 X 10 "nM~ L s- 

0.3s -1 
0.1s -1 

1.00 X 10 -4 nM _1 s" 
0.1s -1 

5.00 X 10 _4 iiM _1 i" 
0.1s -1 

2.00 X 10 -4 rcM -1 s" 

2.00 X 10 -4 rcM -1 s" 

1.00 X 10 -3 rcM -1 s" 
7.00 X 10 -3 n 



1.00 X 10 
2.78 X 10 



-l s -l 

-l s -l 
-l s -l 

-1 -1 



2.84 X 10 -< "nM 
4.41 X 10 -4 rcM -1 s" 

0- 7s -1 

1- 96 X 10 -5 nM -1 s" 
4.8s -1 

1- 06 X 10 -4 nM -1 s" 
2.47 X 10 -3 nM -1 s" 

2- 00 X 10 -4 nM _1 s" 
2-00 X 10 -4 nM -1 s" 



5.70 X 10" 



2.40 X 10" 



The index (p,q,r) in reactions (a) takes values (0,0,0), (1,0,0), (1,0,1), (1,1,0), (2,0,0), (2,0,1), (2,0,2), (2, 1,0), (2, 1,1) 
and (2,2,0). In reactions (b) it takes values (1,0,0), (2,0,0), (2,0,1), (2, 1,0), (3, 0,0), (3,0,1), (3,0, 2), (3, 1,0), (3, 1,1) and 
(3,2,0), while it takes values (2, 2,0), (3, 2,0), (3, 2,1) and (3,3,0) in reactions (c). 
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XIAP Smac 

Apoptosis 



Figure 2: The Fas-induced apoptotic pathway, including two channels. 



in the downstream process. 

According to the law of mass action, the dynamics of the ISS is governed by 28 ordinary differential 
equations 

d § = Q(U). (3.1) 

Here U = U(t) is a column vector with 28 components representing the concentrations of all the 28 
species in the ISS: 

U = {[Casp%l], [Casp%* 2 : CaspZ], [Casp%* 2 : Bid], [Bid], [tBid], [tBid : Bax], [tBid : Bax 2 ], 
[Bcl 2 : tBid], [Bax], [Bcl 2 ■ Bax], [Bcl 2 ], [Cyto.c], [Cyto.c*], [Cyto.c* : Apaf : ATP], 
[Cyto.c* : Apaf : ATP : Casp9], [Cyto.c* : Apaf : ATP : Casp9 2 ], [Apaf], [Casp9% 
[Casp9], [CaspZ], [Casp9* : CaspS], [CaspT], [Smac], [Smac*], [XIAP], 
[Smac* : XIAP], [Casp9 : XIAP], [CaspT : XIAP]) T , 

each element of the vector- valued function Q(U) of U is the change rate of concentration for the corre- 
sponding species 

Q(U) = (-v 7 + v s - v 9 + Via + v ,v 7 - v 8 ,v g - v w ,-v 9 ,v 10 - v n -v 25 ,vn - Vi 2 ,Vi 2 ,V 25 , 

-Vu - V 12 - V 24 , V 2 4, -V 24 - V 25 , ~V 15 , V 15 - V le ,Vi 6 - Vn, Vn - Vis + "19, 
VlS - «19, -Vl6, V19 - V 20 + U 2 l, -Vn - Vis - V 2 2, -V7 - v 20 ,v 20 - V 2 i, 
V& + V21 - V 23 , -l'i 3 ,V 13 - Vu, -V U - V 22 - V 23 ,Vi4,V 22 ,V 23 ) 
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Table 2: Non-zero initial concentrations of the species in the ISS model (Hua. et al. 2005) 



Initial col 



'cntrat: 



Cyto.c 

Smac 

XIAP 

Caap9 

ATP 

Apaf 



Bid 
Bcl2 
Bax 



200.00 
25.00 
75.00 
83.33 
100.00 
100.00 
30.00 
20.00 



10000.00 
100.00 



with Vi(i = 7 ■ ■ ■ 25) the rate of the i-th reaction in Table 1: 



V-i - 


= k 7 [Casp8* 2 ][Casp3] - k- 7 [Casp8* 2 : CaspS], 






v 8 - 


= k$[Casp8* : CaspS], 






v 9 - 


= k 9 [Casp8* 2 ][Bid] - k- g [Casp8* 2 : Bid], 








= kio[Casp8* : Bid], 






I'll 


— ku [tBid] [Bax] — k-n[tBid : Bax], 






Vu 


= ki2[tBid : Bax][Bax] — k-i2[tBid : Bax2], 






V13 


= ki3 [Smac] [tBid : Bax2] 






Vu 


= k u [Smac*][XIAP] - k-u[Smac* : XIAP], 






V15 


— ki^[Cyto.c][tBid : Bax2], 






VlQ 


= k 16 [Cyto.c*][Avaf][ATP] - k m \Cyto.c* : Apaf : 


ATP], 




VVT 


= kn[Cyto.c* : Apaf : ATP][Casp9] - k_n[Cyto.c* 


: Apaf : ATP : 


Casp9], 


Vis, 


= k ls [Cyto.c* : Apaf : ATP : Casp9] [Casp9] - k-is 


[Cyto.c* : Apaf 


: ATP : Casp9 2 


''19 


= ki 9 [Cyto.c* : Apaf : ATP : Casp9 2 ], 






V20 


= k2o[Casp9*][Casp3] — k^2a[Casp9* : CaspS], 






V-21 


= k2i[Casp9* : CaspS], 






V22 


= k 22 [Casp9][XIAP] - k_22[Casp9 : XIAP], 






V23 


= k 23 [Casp3*][XIAP] - k- 2 a[Casp3* : XIAP], 






V24 


= k2i[Bcl2][Bax] — k-2i[Bcl2 ■ Bax], 






V25 


= k25[Bcl 2 ][tBid] - k-25[Bcl 2 : tBid], 







vq is the constant rate of generation for Casp82 from the upstream process and its value was suggested 
in [TH] as vq — O.OOlnMs -1 . In addition, the non-zero initial concentrations for U are taken as in [SI US] 
and are given in Table [H 

In [19] . Okazaki et al. claimed that Smac and XIAP have little effect on the reaction process of the 
ISS by studying the M-D transition behavior of both the ISS model and ISS(wo/S,X) (without Smac and 
XIAP) model. So they did not consider reactions (H13), (H14), (H22) and (H23) in their simplification. 
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However, in our previous paper |23) we found that Smac and XIAP should not be ignored because the 
numerical results of these two models are quite different if initial concentrations are changed. Therefore, 
our sequel discussion will base on the entire ISS system. 

We conclude this section by explaining the M-D transition behavior. It means a D-channel and 
M-channel switching behavior and the quantity of [CaspSj] is a control parameter. When a large amount 
of Casp82 is activated from the upstream process, it will directly induce cell death through the D-channel; 
otherwise, the M-channel plays more important role for cell death. In [B], the authors claimed that the 
effects of D-channel and M-channel can be altered by varying the amount of Casp82 generated by DISC 
|25[ 1261 [25] . which is consistent with previous experiments. 



4 Model reduction 

In this section we use the PEA method to investigate the large apoptosis system p. ID . This system 
contains 13 reversible reactions: (H7), (H9), (Hll), (H12), (H14), (H16), (H17), (H18), (H20), (H22), 
(H23), (H24) and (H25). In our preliminary work [33], we showed that the six reactions (Hll), (H12), 
(H16), (H17), (H24) and (H25) are fast. Because the PEA method has a solid mathematical basis, it 
can be used as a tool to determine whether or not a reversible reaction is fast. After many attempts by 
assuming some of the rest 7 reactions to be fast too, we find that the nine reversible reactions (H7), (H9), 
(Hll), (H12), (H16), (H17), (H18), (H24) and (H25) can be well regarded as fast. 

Here we derive the corresponding simplified model under the assumption that the nine reversible 
reactions are fast. According to the framework in Section 2, we decompose the concentration vector U as 



U = 



Y 

\ Z ) 



Here X stands for the products of the nine reactions, Y for the reactants and Z for the rest: 



X = 



Y = 
Z = 



([Casp8*. ■ Caspi], [Casp8*, '■ Bid], [tBid : Bax], [tBid : Bax2], 
[Bch : Bax], [Bch : tBid], [Cyto.c* : Apaf : ATP], 
[Cyto.c* : Apaf : ATP : Casp% [Cyto.c* : Apaf : ATP : Casp9 2 ]) T ', 
([Casp8*,], [Bid], [tBid], [Bax], [Bch], [Cyto.c*], [Apaf], [Casp9], [Casp3]) T , 
([Cyto.c], [Casp9*], [Casp9* : CaspS], [Casp3% [Smac], [Smac]*, [XIAP], 
[Smac* : XIAP], [Casp9 : XIAP], [CaspT : XIAP]) T . 



(4.1) 



With this decomposition, the kinetic equations in (|3.1[1 can be rewritten as 

^ = -Q 1 (X,Y) + Q 1 (X,Y,Z) 
dt e 

dY 1 a 

— = -Q 2 (X,Y) + Q 2 (X,Y,Z) . 



dt 
dZ 

~dt 



Q 3 (X,Y,Z). 



(4.2) 
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Here the small parameter e characterizes the fastness of the reversible reactions as in Section 2, 

Ql(X, Y) = e(v 7 , Vg, Vu ~ "12, "12, "24, "25, "16 ~ "17, "17 - V 1S ,V 18 ) T , 

Q 2 (X, Y) = e(-v 7 - vg, -vg, -vu - "25, -vu - "12 - V24, -"24 - "25, -vie, -vi6, -vir - vis, ~v 7 ) T , 
stand for the change rates of concentration due to the rapid reactions, and 

Qi(X,Y,Z) = (-v s ,-v 10 ,0,0,0,0,0,v 19 ,-v 19 ) T , 

Q 2 {X, Y, Z) = (v s + vio + v , 0, uio, 0, 0, v i5 , 0, ~v 22 , -v 20 ) T , 

Q 3 (X,Y, Z)= ( - Ui5,«i9 - V 20 + "21, "20 ~ V 2 l,Vg + V 2 1 - "23, -"13, "13 - "14, ~"14 - V 2 2 - "23 , "14, "22 , "23) ' 

are those for the other reactions. It is direct to check that 



Q 2 (X,Y) + eQi(XX) = o 



(4.3) 



with 6 the following constant 9x9-matrix 

/ iiooooooo\ 

010000000 
1 1 1 
1 2 1 
9= 000011000 
1 1 1 
1 1 1 
000000012 
v 100000000y 

Recall that 

v 7 = k 7 [Casp8* 2 ][Casp3] - k- 7 [Casp8* 2 : CaspS], 
v 9 = k 9 [Casp8* 2 ][Bid] - k-g[Casp8* 2 : Bid], 
vu = fcn [tBid] [Bax] — k-u[tBid : Bax], 
V12 = ki 2 [tBid : Bax][Bax] — k-i2\tBid : Bax2], 
"16 = ki 6 [Cyto.c*][Apaf][ATP] - k- W [Cyto.c* : Apaf : ATP], 
vn = k l7 [Cyto.c* : Apaf : ATP][Caps9] - k- 17 [Cyto.c* : Apaf : ATP : Casp9], 
v ls = k ls [Cyto.c* : Apaf : ATP : Casp9] [Casp9] - k- W [Cyto.c* : Apaf : ATP : Casp9 2 ], 
v 2 4 = k24{Bcl 2 ][Bax] — k-24[Bcl 2 ■ Bax], 
"25 = &25 [Bch] [tBid] - k-2b[Bd2 ■ tBid]. 



(4.4) 



Then for any given Y = ([Casp8%\, [Bid], [tBid], [Bax], [Bch], [Cyto.c*], [Apaf], [Casp9], [Caspi]) T with 
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positive components, we use (14.3j) to get 



X = {[Casp8* 2 : Casp3], [Casp8* 2 : Bid], [tBid : Bax], [tBid : Bax 2 ], [Bd 2 ■ Bax], [Bcl 2 ■ tBid], 

[Cyto.c* : Apaf : ATP], [Cyto.c* : Apaf : ATP : Casp% [Cyto.c* : Apaf : ATP : Casp9 2 }) T 

with positive components such that vj = vg — vn — v\ 2 — vie — vyj — uig = «24 = i'25 = 0. Thus, the 
principle of detailed balance is verified. 

After verifying the conditions for the PEA method to be reliable, we turn to write down the simplified 
model. In view of ([3~5]) . we define Y = Y + QX. Then the ODEs in (|3~2"|) become 

^ = -Q 1 (X,Y) + Q 1 (X,Y,Z), 
dt e 

dY 

— = Q 2 (X,Y,Z) + QQ 1 (X,Y,Z), (4.5) 
dt 

^ = Q 3 (X,Y,Z). 

Guided by the framework in Section 2, we solve Qi(X, Y) = 0, namely, 

Vf = k-j[C 'asp8*][C 'asp3] — k- 7 [Casp8* : Caspi] = 0, 
v 9 = k 9 [Casp8* 2 ][Bid] - k-g[Casp8* : Bid] = 0, 
fn = kn[tBid][Bax] — k-n[tBid : Bax] = 0, 
V12 = k\ 2 [tBid : Bax][Bax] — k^\ 2 [tBid : Bax 2 ] = 0, 
< vie = k\Q [Cyto.c*] [Apaf] [ATP] - k^ 16 [Cyto.c* : Apaf : ATP] = 0, 

Vi 7 = k u [Cyto.c* : Apaf : ATP][Caps9] - k- 17 [Cyto.c* : Apaf : ATP : Casp9] = 0, 

v 18 = k ls [Cyto.c* : Apaf : ATP : Casp9] [Casp9] - k^ ls [Cyto.c* : Apaf : ATP : Casp9 2 ] = 0, 

i>24 = k 2 4[Bcl 2 ] [Bax] — k- 2 i[Bcl 2 : Bax] = 0, 

v 25 = k 25 [Bcl 2 ][tBid] - k_ 25 [Bcl 2 : tBid] = 0. 
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From these algebraic equations we can easily solve X in terms of Y 

kj[C asp8*][C asp3] 



[Casp8* 2 : CaspZ] 
[Casp8* 2 : Bid] 
[tBid : Bax] 
[tBid : Bax2] 
[Cyto.c* : Apaf : ATP] 
[Cyto.c* : Apaf : ATP : Casp9] 

[Cyto.c* : Apaf : ATP : Casp%] 

[Bcl 2 ■ Bax] 
[Bcl 2 : tBid] 



fe-7 

k 9 [Casp8*][Bid] 



fell [t Bid] [Bax] 
fc-11 



= K 7 [Casp8*][Casp3], 
-- K 9 [Casp8*][Bid], 
K n [tBid][Bax], 



k la \tBi*Bax][Bax] = Kl2Kn [ tB id][Bax][Bax], 
k 16 [Cyto.c*\[Apaf][ATP\ = Kw { C yto.C*][Apaf][ATP], 

fcir [Cyto.c* -Apaf -ATP] [Casp9] 
fc-17 

Ki 7 K w [Cyto.c*] [Apaf] [ATP] [Casp9] , 

fcis [Cyto.c* : Apaf: AT P:Casp9] [Casp9] 
fc-is 

K 1S K 17 K 16 [Cyto.c*] [Apaf] [ATP] [CaspQ] 2 , 

k 2l [Bcl 2 ][Bax\ = K2 ^ Bd ^[Bax], 

= K 25 [Bcl 2 ][tBid] 



(4.6) 



fc-2 



k 25 [Bcl 2 ][tBid] 

fc-25 



with K t = k, L jk^ L for i = 7, 9, 11, 12, 16, 17, 18, 24, 25. Denote these relations by X = $(Y). 

It is remarkable that the relations rely only on the 9 constants Ki, instead of the 18 constants k±i 
The latter are often not reliably known. 

Substituting X = ^(Y) into the second and third equations in (|4.5I) . we obtain 

dY 

— = Q 2 (*(F), Y, Z) + 0Qi(*(y), Y, Z), 



dZ 
~dt 



= Q 3 (V(Y),Y,Z) 



and thereby gain equations for Y: 



dY dY 

°— = {I + m(Y) Y y l — = (I + m{Y) Y )- l (Q 2 {%Y),Y,Z) + ®Q 1 {*{Y),Y,Z)). 
dt dt 

Consequently, the original system (|3.1j) of 28 ODEs can be approximated by the following 19 ODEs 

dY 



, t = (i + e*(r) y )- 1 (g 2 ($(r), y, z) + oq 1 ^(y),y, z)), 

¥L = Q 3 {V(Y),Y,Z) 



(4.7) 



together with nine algebraic relations 



x = *(y) 



being detailed in (|4.6p . Recall that Y and Z are defined in (|4.1[) . We call this new simplified model as 
ISS-2. 
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Table 3: The ISS skeleton model due to Okazaki et al. (2008) 



Reaction Rate constant 

(SI) CaspS* + Casp3 -CaspS* + CaspS* 6.25.00 X 1 - 6 n M ~ 1 s _ 1 

k a \CaspS%]\Bid] , 

(S2a) CoapS* + Bid -CaspS* + 0.0328tBid : Bai 2 «S2a = [CaspS* ] + K ' fca = i s « = 20nM) 

(S2b) Cyto.c + tBid : Bax 2 -0 . 8G7Cy t o . c* : Apaf : ATP + tBid : Bax 2 1 X 10~ 3 tiM -1 s -1 

(S2c) Cyto.c* : Apaf : ATP + 2Casp9 -Cyto.c* : Apaf : ATP + CaspQ + Casp9* 1.46 X 1 ~ 6 ™ M ~ 1 s - 1 

(S2d) Coap9* + Casp3 -CaspB* + Casp3* 1.96 X 10 — 5 nM~ 1 3~ 1 



5 Numerical simulations 

The purpose of this section is to show the reliability of our ISS-2 model by resorting to numerical 
simulations. Precisely, we compare the ISS-2 model f|4.6[) — f|4. 7[) with the entire ISS model (|3.1j) and 
Okazaki et al.'s ISS skeleton model [H] in several aspects, including the accuracy, M-D transition behavior 
and sensitivity. For the reader's convenience, the skeleton model is given in Table [3] The simulations 
were carried out with Matlab. 



5.1 Accuracy of the ISS-2 model 

We compute the concentration of each species as functions of time t for the entire ISS model, the ISS 



skeleton model and our ISS-2 model, with initial concentrations from Table [5J Fig. 3(a) displays three 
curves of Casp3* as functions of time t corresponding to the three models. In this figure, the equilibrium 
value of Casp3* — the indicator of apoptosis — of the ISS-2 model is almost the same as that of the ISS 
model, whereas that of the skeleton model is slightly larger. The equilibrium values of all other species 
for the ISS model and our ISS-2 model are also very close to each other. The curves of Casp3, Casp9* and 
Bid as functions of time t are given in Fig. |3(b)| Fig. 3(c) and Fig. |3(d)| respectively. These numerical 
results show that our ISS-2 model is a reliable simplification of the entire ISS model. 



5.2 M-D transition behavior 

The M-D transition behavior is explained at the end of Section 3. In [19 , it was reported that initial 
concentrations of Casp9 also have considerable impacts to the transition behavior. In order to study this 
behavior, Okazaki et al. introduced two quantities and Vq in [191 . The former was defined as the ratio 
of the net production of Casp3* via the D-channel to its total production, while the latter represents the 
critical value of vq (the generation rate for CaspS^) corresponding to jd — 0.5. Note that, at 7/3 = 0.5, 
the effect of the M-channel is same as that of the D-channel. 

As previously, we use the initial concentrations from Table [5] and numerically solve the three models 



to obtain three curves of 715 as a function of vq. The results are shown in Fig. 4(a) From this figure, we 
see that the curve given by the ISS-2 model matches that by the ISS model quite well and is obviously 
better than that by the skeleton model. 

The curves of Vq as a function of the initial concentration of Casp9 are shown in Fig. |4(b)| From 
this figure we see that when the initial concentrations of Casp9 are small, the values of Vq for the three 
models are almost the same. However, when the initial concentrations of Casp9 are large, the ISS skeleton 
model behaves quite different from the ISS model. But our ISS-2 model still matches the ISS model very 
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Time[hour] 

(a) The curves of Casp3* as functions of t 



-* — ISS 
-6 — Skeleton 
- - ISS-2 




Tlme[hour] 

(b) The curves of Casp3 as functions of t 




Time[hour] 




2 3 4 5 

Time[hour] 



(c) The curves of Casp9* as functions of t 



(d) The curves of Bid as functions of t 



Figure 3: Comparison of the ISS model (red asterisk), the skeleton model (blue open circle) and our ISS-2 
model (green dot). The initial concentrations are from Tabled (a): the curves of Casp3* as functions of 
t. (b): the curves of Casp3 as functions of t. (c): the curves of Casp9* as functions of t. (d): the curves 
of Bid as functions of t. 
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v /nMs _1 [Casp9J /nM 

(a) M-D transition behaviors due to Casp8?; (b) M-D transition behaviors due to Casp9 

Figure 4: Comparison of the ISS model (red asterisk), the skeleton model (blue open circle) and the ISS-2 
model (green dot) for the M-D transition behavior, (a): M-D transition behavior due to CaspSjjS. (b): 
M-D transition behavior due to Casp9. 

well. All these indicate that our ISS-2 model can well describe the actual M-D transition behavior and 
are much better than the ISS skeleton model. 

5.3 Sensitivity Analysis 

Now we present some results on the sensitivity of our ISS-2 model. Because the half-time — the time for 
Casp3* to attain half of its equilibrium value — is an important quantity to characterize how fast a cell 
will die [6j , we compute this quantity for the three models with initial data changed by one or two orders 
of magnitude higher and lower than the baseline values given in Tabled The numerical results are shown 
in Fig. [5j 

From Fig. [5]we see that, like the full apoptosis model due to Hua et al. [B], our ISS-2 model possesses 
the symmetrical or asymmetrical properties of varying each species to the outcome. The result by the 
ISS-2 model is very similar to that by the ISS model. A bit difference is that the half-time for the ISS-2 
model is a little shorter than that for the ISS model, which is same as for the ISS skeleton model. This 
is expected because the assumption of fast reactions slightly speeds up the whole apoptotic process. 

To evaluate our simplified model, we follow our previous work [23] and introduce the quantity 

equilibrium value of CaspS* for ISS-2 
C3 equilibrium value of CaspS* for ISS 

to examine how different are the equilibrium values of Casp3* for the two models when changing initial 
concentrations of a certain species. When initial concentrations of some species are changed, acy will 
likely change too. For a good simplified model, such a quantity should be close to one. 

We compute ctc3* for changing initial concentrations of each species, including Casp3*, by one or 
two orders of magnitude higher and lower than the baseline values as before. The result is given in Fig. 
[51 This result illustrates that ac3* is insensitive to most of initial concentration changes, except a little 
sensitivity for Bcl2 and Bax. In conclusion, our ISS-2 model well retains the main features of the ISS 
model and therefore can be viewed as a reliable simplification to the original ISS model. 
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Figure 5: Sensitivity analysis of the ISS model, the skeleton model and the ISS-2 model. The overexpres- 
sion or knockdown level of each species is changed one or two orders of magnitude of the baseline values 
while the others are unchanged. 



1.4 



1.2 - 



0.4 



0.2- 








vQ Casp9 Bcl2 Bid Bax Cyto.c Apaf Smac XIAP Casp3 

decrease 1 0OX ^^^^ decrease 1 0X ^^^^ baseline value ^^^^ increase 1 0X increase 1 00X 



Figure 6: The change of acy against initial concentrations. The overexpression or knockdown level 
of each species is changed one or two orders of magnitude of the baseline values while the others are 
unchanged. 
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6 Summary 



In this paper, we develop a general framework of the PEA method together with two conditions, under 
which the method can be justified rigorously. These conditions are the fastness assumption and the 
principle of detailed balance for fast reactions as a whole. Under these conditions, we simplify a general 
system of chemical reactions governed by the law of mass action. This simplification clearly has a solid 
mathematical basis. 

Then we follow the general framework and study the ISS (intracellular-signaling subsystem) model 
as the downstream process of the Fas-signaling pathway model proposed by Hua et al. (2005) for human 
Jurkat T cells. Because the framework has a solid mathematical basis, it can be used as a tool to 
determine whether or not a reaction is relatively fast. After many attempts via numerical tests, we found 
that nine of reactions in the ISS model can be well regarded as fast. 

Knowing that the nine reactions are faster than others, we use the justified PEA method and simplify 
the ISS model to derive a so-called ISS-2 model. It is remarkable that the reversible reactions in apoptosis 
obey the principle of detailed balance naturally. With numerical simulations, we compare the ISS-2 model 
with the ISS model as well as Okazaki's ISS skeleton model in several aspects, including the accuracy, M- 
D transition behavior and sensitivity analysis. All the simulations show that the ISS-2 model is reliable. 
In particular, the new model can very well capture the M-D transition behavior of the ISS model at large 
initial concentrations of Casp9 and therefore improves Okazaki's ISS skeleton model considerably (see 
Fig. [4(b)T>. 

At present, we are trying to simplify the upstream process with the justified PEA method. In the 
future, we will also try to simplify the whole process by correctly combining the PEA method and the 
QSSA method. 
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